Some spatial analysis

R25 Modelers and Story Tellers

1 Introduction

In this tutorial, we learn:

• some spatial analyses.

Resources: Analyzing US Census Data Chapter 7 and R for Data Science.

2 Spatial predicates

How do we find Census tracts that intersect with the Los Angeles metro area?

• Find all Census tracts in CA.

library(tigris)
library(tidyverse)
library(sf)
options(tigris_use_cache = TRUE)

# CRS used: NAD27 / California zone VII
ca_tracts <- tracts(state = "CA", cb = TRUE) %>%
  st_transform(4267) %>%
  print()

Simple feature collection with 9109 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -124.4085 ymin: 32.53438 xmax: -114.1304 ymax: 42.00966
Geodetic CRS:  NAD27
First 10 features:
   STATEFP COUNTYFP TRACTCE             AFFGEOID       GEOID    NAME
1       06      037  294610 1400000US06037294610 06037294610 2946.10
2       06      059  011402 1400000US06059011402 06059011402  114.02
3       06      019  001304 1400000US06019001304 06019001304   13.04
4       06      037  273700 1400000US06037273700 06037273700    2737
5       06      037  207501 1400000US06037207501 06037207501 2075.01
6       06      053  010900 1400000US06053010900 06053010900     109
7       06      059  087403 1400000US06059087403 06059087403  874.03
8       06      037  242100 1400000US06037242100 06037242100    2421
9       06      037  408134 1400000US06037408134 06037408134 4081.34
10      06      001  451704 1400000US06001451704 06001451704 4517.04
               NAMELSAD STUSPS         NAMELSADCO STATE_NAME LSAD   ALAND
1  Census Tract 2946.10     CA Los Angeles County California   CT  823808
2   Census Tract 114.02     CA      Orange County California   CT 1136419
3    Census Tract 13.04     CA      Fresno County California   CT 1307803
4     Census Tract 2737     CA Los Angeles County California   CT  690352
5  Census Tract 2075.01     CA Los Angeles County California   CT  134327
6      Census Tract 109     CA    Monterey County California   CT 3736271
7   Census Tract 874.03     CA      Orange County California   CT  630332
8     Census Tract 2421     CA Los Angeles County California   CT  459615
9  Census Tract 4081.34     CA Los Angeles County California   CT  884018
10 Census Tract 4517.04     CA     Alameda County California   CT 2005363
   AWATER                       geometry
1       0 MULTIPOLYGON (((-118.2624 3...
2       0 MULTIPOLYGON (((-117.9128 3...
3       0 MULTIPOLYGON (((-119.7537 3...
4       0 MULTIPOLYGON (((-118.4583 3...
5       0 MULTIPOLYGON (((-118.255 34...
6   15260 MULTIPOLYGON (((-121.3902 3...
7       0 MULTIPOLYGON (((-117.9145 3...
8       0 MULTIPOLYGON (((-118.2358 3...
9       0 MULTIPOLYGON (((-117.9037 3...
10      0 MULTIPOLYGON (((-121.7975 3...

• Find Los Angeles CBSA.

la_metro <- core_based_statistical_areas(cb = TRUE, year = 2020) %>%
  filter(str_detect(NAME, "Los Angeles")) %>%
  st_transform(4267) %>%
  print()

Simple feature collection with 1 feature and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -118.9437 ymin: 32.80142 xmax: -117.4124 ymax: 34.82331
Geodetic CRS:  NAD27
  CSAFP CBSAFP       AFFGEOID GEOID                               NAME
1   348  31080 310M600US31080 31080 Los Angeles-Long Beach-Anaheim, CA
                                       NAMELSAD LSAD       ALAND     AWATER
1 Los Angeles-Long Beach-Anaheim, CA Metro Area   M1 12566940582 2193795620
                        geometry
1 MULTIPOLYGON (((-118.6035 3...

• Plot:

ggplot() +
  geom_sf(data = ca_tracts, fill = "white", color = "grey") +
  geom_sf(data = la_metro, fill = NA, color = "red") +
  theme_void() +
  labs(title = "Census tracts relative to LA CBSA")

• Spatial subsets that intersect:

ca_tracts_intersect <- ca_tracts %>%
  st_filter(la_metro, .predicate = st_intersects)

ggplot() + 
  geom_sf(data = ca_tracts_intersect, fill = "white", color = "grey") + 
  geom_sf(data = la_metro, fill = NA, color = "red") + 
  theme_void() + 
  labs(title = "Census tracts that intersect or border with the Los Angeles CBSA")

• Predicate st_within:

ca_tracts_within <- ca_tracts %>%
  st_filter(la_metro, .predicate = st_within)

ggplot() + 
  geom_sf(data = ca_tracts_within, fill = "white", color = "grey") + 
  geom_sf(data = la_metro, fill = NA, color = "red") + 
  theme_void() + 
  labs(title = "Census tracts within the Los Angeles CBSA")

2.1 Spatial joins and group-wise spatial analysis

library(tidycensus)
library(tidyverse)
library(sf)

# CRS: NAD83(2011) / Texas Centric Albers Equal Area
tx_cbsa <- get_acs(
  geography = "cbsa",
  variables = "B01003_001", # total population
  year = 2019,
  survey = "acs1",
  geometry = TRUE
) %>%
  filter(str_detect(NAME, "TX")) %>%
  slice_max(estimate, n = 4) %>%
  st_transform(6579)

tx_cbsa

Simple feature collection with 4 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1538292 ymin: 7167200 xmax: 2045975 ymax: 7705975
Projected CRS: NAD83(2011) / Texas Centric Albers Equal Area
  GEOID                                            NAME   variable estimate moe
1 19100      Dallas-Fort Worth-Arlington, TX Metro Area B01003_001  7573136  NA
2 26420 Houston-The Woodlands-Sugar Land, TX Metro Area B01003_001  7066140  NA
3 41700        San Antonio-New Braunfels, TX Metro Area B01003_001  2550960  NA
4 12420     Austin-Round Rock-Georgetown, TX Metro Area B01003_001  2227083  NA
                        geometry
1 MULTIPOLYGON (((1681247 760...
2 MULTIPOLYGON (((2009903 730...
3 MULTIPOLYGON (((1538306 729...
4 MULTIPOLYGON (((1664195 732...

pct_hispanic <- get_acs(
  geography = "tract",
  variables = "DP05_0071P", # percent Hispanic
  state = "TX",
  year = 2019,
  geometry = TRUE
) %>%
  st_transform(6579)

pct_hispanic

Simple feature collection with 5265 features and 5 fields (with 11 geometries empty)
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 873007.1 ymin: 6861658 xmax: 2118112 ymax: 8045409
Projected CRS: NAD83(2011) / Texas Centric Albers Equal Area
First 10 features:
         GEOID                                      NAME   variable estimate
1  48113019204 Census Tract 192.04, Dallas County, Texas DP05_0071P     58.1
2  48377950200 Census Tract 9502, Presidio County, Texas DP05_0071P     95.6
3  48029190601 Census Tract 1906.01, Bexar County, Texas DP05_0071P     86.2
4  48355002301  Census Tract 23.01, Nueces County, Texas DP05_0071P     80.6
5  48441012300    Census Tract 123, Taylor County, Texas DP05_0071P     26.4
6  48051970500 Census Tract 9705, Burleson County, Texas DP05_0071P     20.4
7  48441010400    Census Tract 104, Taylor County, Texas DP05_0071P     64.4
8  48201311900   Census Tract 3119, Harris County, Texas DP05_0071P     84.9
9  48113015500    Census Tract 155, Dallas County, Texas DP05_0071P     56.6
10 48217960500     Census Tract 9605, Hill County, Texas DP05_0071P      9.3
   moe                       geometry
1  5.7 MULTIPOLYGON (((1801563 765...
2  3.7 MULTIPOLYGON (((1060841 729...
3  6.7 MULTIPOLYGON (((1642736 726...
4  3.9 MULTIPOLYGON (((1755212 707...
5  6.3 MULTIPOLYGON (((1522575 758...
6  6.1 MULTIPOLYGON (((1812744 736...
7  7.7 MULTIPOLYGON (((1522727 759...
8  5.8 MULTIPOLYGON (((1950592 729...
9  7.5 MULTIPOLYGON (((1778712 763...
10 4.8 MULTIPOLYGON (((1741591 753...

Spatial join:

hispanic_by_metro <- st_join(
  pct_hispanic,
  tx_cbsa,
  join = st_within,
  suffix = c("_tracts", "_metro"),
  left = FALSE # inner-join
) 

hispanic_by_metro

Simple feature collection with 3189 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1538292 ymin: 7167200 xmax: 2045975 ymax: 7705975
Projected CRS: NAD83(2011) / Texas Centric Albers Equal Area
First 10 features:
   GEOID_tracts                                NAME_tracts variable_tracts
1   48113019204  Census Tract 192.04, Dallas County, Texas      DP05_0071P
3   48029190601  Census Tract 1906.01, Bexar County, Texas      DP05_0071P
8   48201311900    Census Tract 3119, Harris County, Texas      DP05_0071P
9   48113015500     Census Tract 155, Dallas County, Texas      DP05_0071P
11  48439102000   Census Tract 1020, Tarrant County, Texas      DP05_0071P
13  48201450200    Census Tract 4502, Harris County, Texas      DP05_0071P
14  48201450400    Census Tract 4504, Harris County, Texas      DP05_0071P
22  48157670300 Census Tract 6703, Fort Bend County, Texas      DP05_0071P
25  48201554502 Census Tract 5545.02, Harris County, Texas      DP05_0071P
27  48121020503  Census Tract 205.03, Denton County, Texas      DP05_0071P
   estimate_tracts moe_tracts GEOID_metro
1             58.1        5.7       19100
3             86.2        6.7       41700
8             84.9        5.8       26420
9             56.6        7.5       19100
11            15.7        6.4       19100
13            10.7        3.9       26420
14            26.8       13.1       26420
22            27.0        5.8       26420
25            13.4        3.0       26420
27            35.3        8.7       19100
                                        NAME_metro variable_metro
1       Dallas-Fort Worth-Arlington, TX Metro Area     B01003_001
3         San Antonio-New Braunfels, TX Metro Area     B01003_001
8  Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
9       Dallas-Fort Worth-Arlington, TX Metro Area     B01003_001
11      Dallas-Fort Worth-Arlington, TX Metro Area     B01003_001
13 Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
14 Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
22 Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
25 Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
27      Dallas-Fort Worth-Arlington, TX Metro Area     B01003_001
   estimate_metro moe_metro                       geometry
1         7573136        NA MULTIPOLYGON (((1801563 765...
3         2550960        NA MULTIPOLYGON (((1642736 726...
8         7066140        NA MULTIPOLYGON (((1950592 729...
9         7573136        NA MULTIPOLYGON (((1778712 763...
11        7573136        NA MULTIPOLYGON (((1746016 762...
13        7066140        NA MULTIPOLYGON (((1925521 730...
14        7066140        NA MULTIPOLYGON (((1922292 730...
22        7066140        NA MULTIPOLYGON (((1935603 728...
25        7066140        NA MULTIPOLYGON (((1918552 732...
27        7573136        NA MULTIPOLYGON (((1766859 768...

hispanic_by_metro %>%
  mutate(NAME_metro = str_replace(NAME_metro, ", TX Metro Area", "")) %>%
  ggplot() + 
  geom_density(aes(x = estimate_tracts), color = "navy", fill = "navy", 
               alpha = 0.4) + 
  theme_minimal() + 
  facet_wrap(~NAME_metro) + 
  labs(title = "Distribution of Hispanic/Latino population by Census tract",
       subtitle = "Largest metropolitan areas in Texas",
       y = "Kernel density estimate",
       x = "Percent Hispanic/Latino in Census tract")

Warning: Removed 9 rows containing non-finite values (`stat_density()`).

hispanic_by_metro

Simple feature collection with 3189 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1538292 ymin: 7167200 xmax: 2045975 ymax: 7705975
Projected CRS: NAD83(2011) / Texas Centric Albers Equal Area
First 10 features:
   GEOID_tracts                                NAME_tracts variable_tracts
1   48113019204  Census Tract 192.04, Dallas County, Texas      DP05_0071P
3   48029190601  Census Tract 1906.01, Bexar County, Texas      DP05_0071P
8   48201311900    Census Tract 3119, Harris County, Texas      DP05_0071P
9   48113015500     Census Tract 155, Dallas County, Texas      DP05_0071P
11  48439102000   Census Tract 1020, Tarrant County, Texas      DP05_0071P
13  48201450200    Census Tract 4502, Harris County, Texas      DP05_0071P
14  48201450400    Census Tract 4504, Harris County, Texas      DP05_0071P
22  48157670300 Census Tract 6703, Fort Bend County, Texas      DP05_0071P
25  48201554502 Census Tract 5545.02, Harris County, Texas      DP05_0071P
27  48121020503  Census Tract 205.03, Denton County, Texas      DP05_0071P
   estimate_tracts moe_tracts GEOID_metro
1             58.1        5.7       19100
3             86.2        6.7       41700
8             84.9        5.8       26420
9             56.6        7.5       19100
11            15.7        6.4       19100
13            10.7        3.9       26420
14            26.8       13.1       26420
22            27.0        5.8       26420
25            13.4        3.0       26420
27            35.3        8.7       19100
                                        NAME_metro variable_metro
1       Dallas-Fort Worth-Arlington, TX Metro Area     B01003_001
3         San Antonio-New Braunfels, TX Metro Area     B01003_001
8  Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
9       Dallas-Fort Worth-Arlington, TX Metro Area     B01003_001
11      Dallas-Fort Worth-Arlington, TX Metro Area     B01003_001
13 Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
14 Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
22 Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
25 Houston-The Woodlands-Sugar Land, TX Metro Area     B01003_001
27      Dallas-Fort Worth-Arlington, TX Metro Area     B01003_001
   estimate_metro moe_metro                       geometry
1         7573136        NA MULTIPOLYGON (((1801563 765...
3         2550960        NA MULTIPOLYGON (((1642736 726...
8         7066140        NA MULTIPOLYGON (((1950592 729...
9         7573136        NA MULTIPOLYGON (((1778712 763...
11        7573136        NA MULTIPOLYGON (((1746016 762...
13        7066140        NA MULTIPOLYGON (((1925521 730...
14        7066140        NA MULTIPOLYGON (((1922292 730...
22        7066140        NA MULTIPOLYGON (((1935603 728...
25        7066140        NA MULTIPOLYGON (((1918552 732...
27        7573136        NA MULTIPOLYGON (((1766859 768...

median_by_metro <- hispanic_by_metro %>%
  group_by(NAME_metro) %>%
  summarize(median_hispanic = median(estimate_tracts, na.rm = TRUE))

median_by_metro

Simple feature collection with 4 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 1538292 ymin: 7167200 xmax: 2045975 ymax: 7705975
Projected CRS: NAD83(2011) / Texas Centric Albers Equal Area
# A tibble: 4 × 3
  NAME_metro                                   media…¹                  geometry
  <chr>                                          <dbl>            <GEOMETRY [m]>
1 Austin-Round Rock-Georgetown, TX Metro Area     25.9 POLYGON ((1700741 730245…
2 Dallas-Fort Worth-Arlington, TX Metro Area      22.6 POLYGON ((1737530 756507…
3 Houston-The Woodlands-Sugar Land, TX Metro …    32.4 MULTIPOLYGON (((1901162 …
4 San Antonio-New Braunfels, TX Metro Area        53.5 POLYGON ((1619499 717051…
# … with abbreviated variable name ¹​median_hispanic

3 Distance and proximity analysis

3.1 Distance

All Level I trauma centers in LA metro area:

library(tigris)
library(sf)
library(tidyverse)
options(tigris_use_cache = TRUE)

# CRS: NAD83 / Iowa North
oc_subdivs <- county_subdivisions(
  state = "CA", 
  county = c("Orange"),
  cb = TRUE, 
  year = 2019
  ) %>%
  st_transform(4267)

hospital_url <- "https://services1.arcgis.com/Hp6G80Pky0om7QvQ/arcgis/rest/services/Hospital/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"

trauma <- st_read(hospital_url) %>%
  filter(str_detect(TRAUMA, "LEVEL I")) %>%
  st_transform(4267) %>%
  distinct(ID, .keep_all = TRUE)

Reading layer `OGRGeoJSON' from data source 
  `https://services1.arcgis.com/Hp6G80Pky0om7QvQ/arcgis/rest/services/Hospital/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson' 
  using driver `GeoJSON'
Simple feature collection with 8013 features and 32 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -176.6403 ymin: -14.29024 xmax: 145.7245 ymax: 71.29773
Geodetic CRS:  WGS 84

oc_trauma <- trauma %>%
  st_filter(
    oc_subdivs, 
    .predicate = st_is_within_distance,
    dist = 10000 # 100KM
    ) %>%
  print()

Simple feature collection with 8 features and 32 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -118.186 ymin: 33.56038 xmax: -117.6645 ymax: 33.80818
Geodetic CRS:  NAD27
  OBJECTID         ID
1      478 0070892868
2      495 0030492691
3      502 0051292668
4      503 0052792705
5      519 0006092691
6      881 0026390806
7      886 0046390801
8      948 0011090806
                                                            NAME
1                           CHILDREN'S HOSPITAL OF ORANGE COUNTY
2                                    PROVIDENCE MISSION HOSPITAL
3                 UNIVERSITY OF CALIFORNIA IRVINE MEDICAL CENTER
4                            ORANGE COUNTY GLOBAL MEDICAL CENTER
5                                 CHILDREN'S HOSPITAL AT MISSION
6                         MEMORIALCARE LONG BEACH MEDICAL CENTER
7                                        ST. MARY MEDICAL CENTER
8 MEMORIALCARE MILLER CHILDREN'S AND WOMEN'S HOSPITAL LONG BEACH
                  ADDRESS          CITY STATE   ZIP ZIP4      TELEPHONE
1      1201 W LA VETA AVE        ORANGE    CA 92868 4203 (818) 267-7478
2 27700 MEDICAL CENTER RD MISSION VIEJO    CA 92691 6426 (949) 364-4840
3       101 THE CITY DR S        ORANGE    CA 92868 3201 (714) 453-7880
4       1001 N TUSTIN AVE     SANTA ANA    CA 92705 3502 (714) 953-3610
5 27700 MEDICAL CENTER RD MISSION VIEJO    CA 92691 6426 (818) 267-7478
6       2801 ATLANTIC AVE    LONG BEACH    CA 90806 1701 (562) 933-1902
7         1050 LINDEN AVE    LONG BEACH    CA 90813 3321 (562) 491-9000
8       2801 ATLANTIC AVE    LONG BEACH    CA 90806 1701 (562) 933-8002
                TYPE STATUS POPULATION      COUNTY COUNTYFIPS COUNTRY LATITUDE
1           CHILDREN   OPEN        334      ORANGE      06059     USA 33.78081
2 GENERAL ACUTE CARE   OPEN        345      ORANGE      06059     USA 33.56102
3 GENERAL ACUTE CARE   OPEN        459      ORANGE      06059     USA 33.78879
4 GENERAL ACUTE CARE   OPEN        254      ORANGE      06059     USA 33.75438
5           CHILDREN   OPEN         54      ORANGE      06059     USA 33.56040
6 GENERAL ACUTE CARE   OPEN        411 LOS ANGELES      06037     USA 33.80820
7 GENERAL ACUTE CARE   OPEN        360 LOS ANGELES      06037     USA 33.78042
8           CHILDREN   OPEN        357 LOS ANGELES      06037     USA 33.80817
  LONGITUDE NAICS_CODE                             NAICS_DESC
1 -117.8658     622110          CHILDREN'S HOSPITALS, GENERAL
2 -117.6654     622110 GENERAL MEDICAL AND SURGICAL HOSPITALS
3 -117.8887     622110 GENERAL MEDICAL AND SURGICAL HOSPITALS
4 -117.8334     622110 GENERAL MEDICAL AND SURGICAL HOSPITALS
5 -117.6658     622110          CHILDREN'S HOSPITALS, GENERAL
6 -118.1869     622110 GENERAL MEDICAL AND SURGICAL HOSPITALS
7 -118.1861     622110 GENERAL MEDICAL AND SURGICAL HOSPITALS
8 -118.1855     622110          CHILDREN'S HOSPITALS, GENERAL
                                                                                                                                                           SOURCE
1 https://data.chhs.ca.gov/dataset/3b5b80e8-6b8d-4715-b3c0-2699af6e72e5/resource/098bbc36-044d-441f-9442-1f4db4d8aaa0/download/healthcare_facility_locations.xlsx
2 https://data.chhs.ca.gov/dataset/3b5b80e8-6b8d-4715-b3c0-2699af6e72e5/resource/098bbc36-044d-441f-9442-1f4db4d8aaa0/download/healthcare_facility_locations.xlsx
3 https://data.chhs.ca.gov/dataset/3b5b80e8-6b8d-4715-b3c0-2699af6e72e5/resource/098bbc36-044d-441f-9442-1f4db4d8aaa0/download/healthcare_facility_locations.xlsx
4 https://data.chhs.ca.gov/dataset/3b5b80e8-6b8d-4715-b3c0-2699af6e72e5/resource/098bbc36-044d-441f-9442-1f4db4d8aaa0/download/healthcare_facility_locations.xlsx
5 https://data.chhs.ca.gov/dataset/3b5b80e8-6b8d-4715-b3c0-2699af6e72e5/resource/098bbc36-044d-441f-9442-1f4db4d8aaa0/download/healthcare_facility_locations.xlsx
6 https://data.chhs.ca.gov/dataset/3b5b80e8-6b8d-4715-b3c0-2699af6e72e5/resource/098bbc36-044d-441f-9442-1f4db4d8aaa0/download/healthcare_facility_locations.xlsx
7 https://data.chhs.ca.gov/dataset/3b5b80e8-6b8d-4715-b3c0-2699af6e72e5/resource/098bbc36-044d-441f-9442-1f4db4d8aaa0/download/healthcare_facility_locations.xlsx
8 https://data.chhs.ca.gov/dataset/3b5b80e8-6b8d-4715-b3c0-2699af6e72e5/resource/098bbc36-044d-441f-9442-1f4db4d8aaa0/download/healthcare_facility_locations.xlsx
    SOURCEDATE    VAL_METHOD     VAL_DATE       WEBSITE  STATE_ID
1 1.633046e+12 IMAGERY/OTHER 1.421107e+12 NOT AVAILABLE 060000007
2 1.633046e+12 IMAGERY/OTHER 1.646006e+12 NOT AVAILABLE 060000060
3 1.633046e+12 IMAGERY/OTHER 1.391990e+12 NOT AVAILABLE 060000071
4 1.633046e+12 IMAGERY/OTHER 1.391990e+12 NOT AVAILABLE 060000073
5 1.633046e+12 IMAGERY/OTHER 1.394582e+12 NOT AVAILABLE 060001207
6 1.633046e+12 IMAGERY/OTHER 1.391990e+12 NOT AVAILABLE 930000095
7 1.633046e+12       IMAGERY 1.391990e+12 NOT AVAILABLE 930000103
8 1.633046e+12 IMAGERY/OTHER 1.394582e+12 NOT AVAILABLE 930001709
                                   ALT_NAME ST_FIPS              OWNER
1      CHILDREN'S HOSPITAL OF ORANGE COUNTY      06         NON-PROFIT
2  MISSION HOSPITAL REGIONAL MEDICAL CENTER      06         NON-PROFIT
3 UNIVERSITY OF CALIFORNIA BOARD OF REGENTS      06 GOVERNMENT - STATE
4 ORANGE COUNTY GLOBAL MEDICAL CENTER, INC.      06        PROPRIETARY
5            CHILDREN'S HOSPITAL AT MISSION      06         NON-PROFIT
6        LONG BEACH MEMORIAL MEDICAL CENTER      06         NON-PROFIT
7                            DIGNITY HEALTH      06         NON-PROFIT
8        LONG BEACH MEMORIAL MEDICAL CENTER      06         NON-PROFIT
  TTL_STAFF BEDS                             TRAUMA HELIPAD
1      -999  334                 LEVEL II PEDIATRIC       Y
2      -999  345 LEVEL II ADULT, LEVEL II PEDIATRIC       Y
3      -999  459  LEVEL I ADULT, LEVEL II PEDIATRIC       Y
4      -999  254                           LEVEL II       Y
5      -999   54                 LEVEL II PEDIATRIC       Y
6      -999  411 LEVEL II ADULT, LEVEL II PEDIATRIC       Y
7      -999  360                           LEVEL II       Y
8      -999  357 LEVEL II ADULT, LEVEL II PEDIATRIC       Y
                    geometry
1 POINT (-117.8649 33.78079)
2 POINT (-117.6645 33.56099)
3 POINT (-117.8878 33.78877)
4 POINT (-117.8325 33.75437)
5  POINT (-117.665 33.56038)
6  POINT (-118.186 33.80818)
7 POINT (-118.1852 33.78041)
8 POINT (-118.1846 33.80816)

ggplot() + 
  geom_sf(data = oc_subdivs, color = "NA", fill = "grey50") + 
  geom_sf(data = oc_trauma, color = "red") + 
  theme_void()

Calculate distances:

dist <- oc_subdivs %>%
  st_centroid() %>%
  st_distance(oc_trauma) 

Warning: st_centroid assumes attributes are constant over geometries

dist[1:5, 1:5]

Units: [m]
          [,1]      [,2]     [,3]      [,4]      [,5]
[1,]  5319.888 34710.384  5175.15  8559.115 34748.184
[2,] 16062.989 16177.958 18312.45 11951.229 16229.072
[3,] 30411.212  5009.441 32216.58 26446.777  4941.636
[4,] 14796.804 36155.001 13460.93 16553.538 36153.336
[5,] 30770.224 17573.176 33060.38 27027.294 17651.429

Histogram of distances to the nearest hospital:

min_dist <- dist %>%
  apply(1, min) %>%
  as.vector() %>%
  magrittr::divide_by(1000)

hist(min_dist)

3.2 Travel time

library(mapboxapi)

# subset_tracts <- slice_sample(lametro_tracts, n = 15)

times <- mb_matrix(
  origins = oc_subdivs, 
  destinations = oc_trauma,
  profile = "driving"
  )

times

         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
[1,] 11.09333 37.41833 13.44333 14.19500 38.25000 32.38167 39.15333 30.74167
[2,] 19.35333 22.70500 21.11833 15.65167 23.53667 38.39500 44.28000 36.75500
[3,] 34.22167 17.82333 35.62333 30.15667 18.65500 45.61667 51.50167 43.97667
[4,] 22.75500 36.13333 22.94833 21.59500 36.96500 22.63667 28.52167 20.99667
[5,] 40.49333 31.63500 42.25833 36.79167 32.46667 59.53500 65.42000 57.89500
[6,] 19.70667 26.48667 21.75167 16.07000 27.31833 31.10167 36.98667 29.46167
[7,] 35.59500 19.23500 37.36000 31.89333 20.06667 54.63667 60.52167 52.99667

min_time <- apply(times, 1, min)

oc_subdivs$time <- min_time

ggplot(oc_subdivs, aes(fill = time)) + 
  geom_sf(color = NA) + 
  scale_fill_viridis_c(option = "magma") + 
  theme_void() + 
  labs(fill = "Time (minutes)",
       title = "Travel time to nearest Level I trauma hospital",
       subtitle = "Census tracts in Orange County, CA",
       caption = "Data sources: US Census Bureau, US DHS, Mapbox")

3.3 Catchment areas: buffers and isochrones

ucla_hosp <- filter(trauma, ID == "0038190095")

# 5 km buffer
buf5km <- st_buffer(ucla_hosp, dist = 5000) 

# 10 min isochrone given current traffic
iso10min <- mb_isochrone(
  ucla_hosp, 
  time = 10, 
  profile = "driving-traffic",
  # depart_at = "2023-07-05T17:00"
  )

Mapping:

library(leaflet)
library(leafsync)

hospital_icon <- makeAwesomeIcon(icon = "ios-medical", 
                                 markerColor = "red",
                                 library = "ion")

# The Leaflet package requires data be in CRS 4326
map1 <- leaflet() %>% 
  addTiles() %>%
  addPolygons(data = st_transform(buf5km, 4326)) %>% 
  addAwesomeMarkers(data = st_transform(ucla_hosp, 4326),
                    icon = hospital_icon)

map2 <- leaflet() %>% 
  addTiles() %>%
  addPolygons(data = iso10min) %>% 
  addAwesomeMarkers(data = st_transform(ucla_hosp, 4326),
                    icon = hospital_icon)

sync(map1, map2)

4 Spatial neighborhoods and spatial weights matrices

library(tidycensus)
library(tidyverse)
library(tigris)
library(sf)
library(spdep)
options(tigris_use_cache = TRUE)

# CRS: 4267
la_cbsa <- core_based_statistical_areas(cb = TRUE, year = 2020) %>%
  filter(str_detect(NAME, "Los Angeles")) %>%
  st_transform(4267)

la_tracts <- get_acs(
  geography = "tract",
  variables = "B01002_001",
  state = "CA",
  year = 2020,
  geometry = TRUE
) %>%
  st_transform(4267) %>%
  st_filter(la_cbsa, .predicate = st_within) %>%
  na.omit()

ggplot(la_tracts) + 
  geom_sf(aes(fill = estimate), color = NA) + 
  scale_fill_viridis_c() + 
  theme_void()

4.1 Spatial neighborhoods

neighbors <- poly2nb(la_tracts, queen = TRUE)

summary(neighbors)

Neighbour list object:
Number of regions: 3088 
Number of nonzero links: 19730 
Percentage nonzero weights: 0.2069057 
Average number of links: 6.389249 
Link number distribution:

  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  18  21  26  38 
  4  12  70 266 606 766 693 368 174  76  28  12   3   4   1   1   1   1   1   1 
4 least connected regions:
1076 1347 1815 2356 with 1 link
1 most connected region:
798 with 38 links

la_coords <- la_tracts %>%
  st_centroid() %>%
  st_coordinates()

Warning: st_centroid assumes attributes are constant over geometries

plot(la_tracts$geometry)
plot(neighbors, 
     coords = la_coords, 
     add = TRUE, 
     col = "blue", 
     points = FALSE)

Spatial weight matrix:

weights <- nb2listw(neighbors, style = "W")

weights$weights[[1]]

[1] 0.25 0.25 0.25 0.25

5 Global and local spatial autocorrelation

First law of geography > Everything is related to everything else, but near things are more related than distant things.

la_tracts$lag_estimate <- lag.listw(weights, la_tracts$estimate)
la_tracts

Simple feature collection with 3088 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -118.9437 ymin: 32.80142 xmax: -117.4124 ymax: 34.82331
Geodetic CRS:  NAD27
First 10 features:
         GEOID                                                 NAME   variable
1  06037300100    Census Tract 3001, Los Angeles County, California B01002_001
2  06037620601 Census Tract 6206.01, Los Angeles County, California B01002_001
3  06037291210 Census Tract 2912.10, Los Angeles County, California B01002_001
4  06037573902 Census Tract 5739.02, Los Angeles County, California B01002_001
5  06037571502 Census Tract 5715.02, Los Angeles County, California B01002_001
6  06037570003 Census Tract 5700.03, Los Angeles County, California B01002_001
7  06037603301 Census Tract 6033.01, Los Angeles County, California B01002_001
8  06037603500    Census Tract 6035, Los Angeles County, California B01002_001
9  06037554516 Census Tract 5545.16, Los Angeles County, California B01002_001
10 06037554517 Census Tract 5545.17, Los Angeles County, California B01002_001
   estimate moe                       geometry lag_estimate
1      44.5 3.9 MULTIPOLYGON (((-118.2477 3...     45.15000
2      38.7 3.3 MULTIPOLYGON (((-118.3718 3...     39.00000
3      37.9 3.3 MULTIPOLYGON (((-118.2909 3...     39.48571
4      51.4 3.2 MULTIPOLYGON (((-118.0807 3...     41.05000
5      42.4 3.6 MULTIPOLYGON (((-118.1918 3...     35.37778
6      39.6 3.1 MULTIPOLYGON (((-118.1458 3...     35.64286
7      45.5 5.0 MULTIPOLYGON (((-118.3126 3...     44.06667
8      40.2 2.1 MULTIPOLYGON (((-118.3256 3...     38.05000
9      47.2 5.4 MULTIPOLYGON (((-118.0714 3...     43.44444
10     39.1 3.7 MULTIPOLYGON (((-118.0628 3...     46.78333

ggplot(la_tracts, aes(x = estimate, y = lag_estimate)) + 
  geom_point(alpha = 0.3) + 
  geom_abline(color = "red") + 
  theme_minimal() + 
  labs(title = "Median age by Census tract, Los Angeles CA",
       x = "Median age",
       y = "Spatial lag, median age", 
       caption = "Data source: 2016-2020 ACS.\nSpatial relationships based on queens-case polygon contiguity.")

Moran’s I test:

moran.test(la_tracts$estimate, weights)

    Moran I test under randomisation

data:  la_tracts$estimate  
weights: weights    

Moran I statistic standard deviate = 46.001, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.4723661752     -0.0003239391      0.0001055896 

5.1 Local spatial autocorrelation

# For Gi*, re-compute the weights with `include.self()`
localg_weights <- nb2listw(include.self(neighbors))

la_tracts$localG <- localG(la_tracts$estimate, localg_weights)

ggplot(la_tracts) + 
  geom_sf(aes(fill = localG), color = NA) + 
  scale_fill_distiller(palette = "RdYlBu") + 
  theme_void() + 
  labs(fill = "Local Gi* statistic")

la_tracts <- la_tracts %>%
  mutate(hotspot = case_when(
    localG >= 2.576 ~ "High cluster",
    localG <= -2.576 ~ "Low cluster",
    TRUE ~ "Not significant"
  ))

ggplot(la_tracts) + 
  geom_sf(aes(fill = hotspot), color = "grey90", size = 0.1) + 
  scale_fill_manual(values = c("red", "blue", "grey")) + 
  theme_void()

5.2 Identifying clusters and spatial outliers with local indicators of spatial association (LISA)

set.seed(1983)

la_tracts$scaled_estimate <- as.numeric(scale(la_tracts$estimate))

la_lisa <- localmoran_perm(
  la_tracts$scaled_estimate, 
  weights, 
  nsim = 999L, 
  alternative = "two.sided"
) %>%
  as_tibble() %>%
  set_names(c("local_i", "exp_i", "var_i", "z_i", "p_i",
              "p_i_sim", "pi_sim_folded", "skewness", "kurtosis"))

la_lisa_df <- la_tracts %>%
  select(GEOID, scaled_estimate) %>%
  mutate(lagged_estimate = lag.listw(weights, scaled_estimate)) %>%
  bind_cols(la_lisa)

Cluster:

la_lisa_clusters <- la_lisa_df %>%
  mutate(lisa_cluster = case_when(
    p_i >= 0.05 ~ "Not significant",
    scaled_estimate > 0 & local_i > 0 ~ "High-high",
    scaled_estimate > 0 & local_i < 0 ~ "High-low",
    scaled_estimate < 0 & local_i > 0 ~ "Low-low",
    scaled_estimate < 0 & local_i < 0 ~ "Low-high"
  ))

LISA quadrant plot:

color_values <- c(`High-high` = "red", 
                  `High-low` = "pink", 
                  `Low-low` = "blue", 
                  `Low-high` = "lightblue", 
                  `Not significant` = "white")

ggplot(la_lisa_clusters, aes(x = scaled_estimate, 
                              y = lagged_estimate,
                              fill = lisa_cluster)) + 
  geom_point(color = "black", shape = 21, size = 2) + 
  theme_minimal() + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed") + 
  scale_fill_manual(values = color_values) + 
  labs(x = "Median age (z-score)",
       y = "Spatial lag of median age (z-score)",
       fill = "Cluster type")

ggplot(la_lisa_clusters, aes(fill = lisa_cluster)) + 
  geom_sf(size = 0.1) + 
  theme_void() + 
  scale_fill_manual(values = color_values) + 
  labs(fill = "Cluster type")